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Abstract 

We present an implementation of realistic static friction in molecular dynamics (MD) 
simulations of granular particles. In our model, to break contacts between two particles, 
one has to apply a finite amount of force, determined by the Coulomb criterion. Using a two 
dimensional model, we show that piles generated by avalanches have a finite angle of repose 
9r (finite slopes). Furthermore, these piles are stable under tilting by an angle smaller 
than a non-zero tilting angle 9t, showing that On is different from the angle of marginal 
stability 0ms ^ which is the maximum angle of stable piles. These measured angles are 
compared to a theoretical approximation. We also measure 9ms by continuously adding 
particles on the top of a stable pile. 

PACS numbers: 46.10+z, 62.20-x 
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1. Introduction 



Systems of granular particles (e.g. sands) exhibit many interesting phenomena. 
The formation of spontaneous heap^~^ and convection cells under vibration, density 
waves found in the outflow through hoppers^^"^^ and segregation of particles^^~^° are just 
a few examples. These phenomena are consequences of the unusual dynamical response of 
the system. One of the characteristic property of granular systems is that it can behave 
like both a solid and a fluid. One can pour (like a fluid) sand grains on a table, and they 
form a stable pile with finite slope (like a solid). Part of the reason why it acts like a 
solid is due to static friction. By static friction, we mean that one has to apply a force 
larger than certain threshold in order to break a contact between particles. This threshold 
is determined by, for example, the Coulomb criterion. Static friction is responsible for 
many static structures (e.g. sand pile), and have a possible implication in the dynamics of 
granular systems. Despite its importance, the effect of static friction has been much less 
studied as compared to other microscopic mechanisms. This is mainly due to the difficulty 
of including static friction to a theoretical framework or a simulational scheme. 

In this paper, we present an implementation of static friction in a molecular dynamics 
(MD) simulation, which uses a scheme introduced by Cundall and Strack.^^ Using this 
code, we generate piles by first filling a (two dimensional) box with grains, then removing 
a sidewalk The slope of the pile is finite, which is related to the finite "angle of repose 
(Or)." Here, tanOR is defined to be the slope of the pile. This angle is strongly dependent 
on the friction coefficient //, and rather insensitive to other parameters of the system. 
Furthermore, we find that the pile obtained above is stable under tilting by an angle 
smaller than the finite tilting angle 9t, where 9t is typically a few degrees. This suggests 
that the angle of marginal stability OmSj the maximum angle of stable piles, is larger 
than Or, which has been observed for real sandpile experiments. ^'^^ We also study the 
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situation of keeping on adding particles on a stable pile. The angle, at which the pile 
becomes unstable, can be interpreted as 9ms- We also propose a theoretical method to 
calculate 6ms and 9t from an approximate stress distribution obtained by Liffman et al?^ 
The theoretical values show similar dependency of the measured angles on the friction 
coefficient ji. 

2. Definition of the model 

The interaction between real sand grains is too complicated to construct a model, by 
which all the properties of a granular system are described accurately. Instead of construct- 
ing a model to reproduce all the details from the beginning, it is often advantageous if one 
identifies the basic ingredients of the system, construct a model with these ingredients. It 
is often true that the qualitative behavior of a system is independent of the fine details of 
the model. Some important ingredients for a granular system are (1) repulsion between 
two particles in contact, (2) dissipation of energy during collision. In certain cases, the 
rotation of the particles could be important. In the previous molecular dynamics (MD) 
simulations of granular systems, most of these ingredients, if not all, are incorporated. For 
example, the repulsion and dissipation is included in most MD simulations of granular 
particles. A few of these simulations also included the rotational degree 
of freedom. -^^'-^^'^^'^^"^^ Here, we will construct a model which includes the repulsion, dis- 
sipation and static friction. But, the model does not have rotational degrees of freedom. 

An individual grain is modeled by a spherical particle. These particles interact with 
each other only if they are in contact. Consider two particles i and j in contact in two 
dimenstion. Let the coordinate for the center of particle i (j) be Ri (Rj), and f= Ri — Rj. 
A vector n is defined to be a unit vector parallel to the line joining the centers of two 
particles, f/r. Another vector s, which is orthogonal to n, is obtained by rotating n by 
7r/2 in clockwise direction. We also define the relative velocity v to be Vi — Vj, and the 
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radius of particle i to be a^. 

The force on particle i exerted by particle j, Fj-,i, can be written as 

Fj^i=Fnfi+FsS, (1) 

where the normal force F„ is given by 

Fn = kn{ai + ttj - rfl"^ - 'jn'meiv-n). (2a) 

The first term of (2a) is the three-dimensional Hertzian repulsion due to the elastic de- 
formation, where is the elastic constant of the material. The second term is a ve- 
locity dependent friction term, which is introduced to dissipate energy from the system. 
Here, 7^ is a constant controlling the amount of dissipation, and mg is the effective mass 
mirrij/ {rrii -\- rrij). The second term of (1), the shear force Fg is 

Fg = -7sme(v • s) - sign(5s) • min(/cs(5s, /iFn), (26) 

where the first term is a velocity dependent damping term similar to the one in (2a). 
The second term of (2b) simulates the static friction. The basic idea is the following. 
When two particles start to touch each other, one puts a "virtual" spring in the shear 
direction. For the total shear displacement 5s during the contact, there is a restoring 
force, kgSs, which is a counter-acting frictional force. The maximum value of this restoring 
force is given by Coulomb's criterion — ^F^. When particles are no longer in contact with 
each other, we remove the spring. We want to emphasize that one has to know the total 
shear displacement of particles during the contact, not the instantaneous displacement to 
calculate the static friction. In other words, one has to remember whether a contact is 
new or old. The rotation of the particles is not included in the present simulation. 

The particles can also interact with walls. If particle i is in contact with a wall, the 
force exerted by the wall on the particle has exactly the same form as eqs. (2) with aj = 
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and rrij = oo. Also, there is a gravitational field. The force acting on particle i by the 
field is —rriig. The total force acting on particle i is the vector sum of the particle-particle 
interaction(s), the particle- wall interaction(s) and the gravitational force. 

The trajectory of a particle is calculated by the fifth order predictor-corrector method.^* 
We use two Verlet tables. One is a usual table with finite skin. The other table is a list 
of pairs of actually interacting particles, which is used to calculate the static friction term. 
For a typical situation, the CPU time needed to run 1372 particles is about 0.01 seconds 
per iteration on a Cray-YMP, which is comparable to the speed of the layered-link-cell 
implementation of a short range Lennard-Jones system.^^'^° 



3. Obtaining the angle of repose 

As a non-trivial check whether the static friction term is working, we measure Or as 
follows. We start by randomly putting particles in a box of width W and height H. The 
particles fall down due to gravity, and loose their energy due to dissipation. After a long 
time, they fill the box with no significant motion. We show, in Fig. 1(a), an example of the 
system at this stage. Here, parameters are /x = 0.2, kn = 10^, kg = 10"^, 7^ = 5 x 10^, 7s = 
7„/100, and for walls, k^, is chosen to be 2 x 10^. We checked the motions of the particles by 
monitoring the total kinetic energy of the system. The average kinetic energy per particle 
during the whole sequence is shown in Fig. 1(b). The kinetic energy sharply rises when the 
wall is removed. The pile relaxes in an oscillatory manner (see Fig. 1(b)). Next, we remove 
the right wall, let the particles move out of the box, and wait until the system reaches a 
new equilibrium. The figure 1(c) is the equilibrium reached by starting from Fig. 1(a). As 
shown in the figure, the new state has non-zero slope. 

We try a few ways to measure 9r. We first divide the box into several vertical cells, 
the width of each cell is equal to the average diameter of the particles. For the center of 
each cell we find the maximum position h{x) of the particles in the cell. The line joined 
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by these positions is a "surface" of the structure. Having determined h{x), we use three 
different ways to measure the slope: (1) By joining h{0) and h{W), (2) by fitting a straight 
line to h{x) using the method of least squares, and (3) by fitting a parabola to Yli=o ^(^) 
by the least squares method. Here, if h(i) is a straight line, the sum is a parabola. In the 
case of h{x) being a straight line, these three methods should give identical results. In 
our simulation, the slopes obtained by different methods differ from each other by a few 
degrees. For example, we obtain (1) 20.14 ± 2.15, (2) 18.90 ± 1.49 and (3) 17.88 ± 1.76 
for 400 particle system with /j, = 0.2. Here the angles are averaged over 10 samples. We 
find, on the average, the angle by method (1) is larger than (2), and (2) is larger than 
(3), although they are within the error bars of each other. It is quite possible that these 
differences come from the finite size of the system. Prom now on, we use only the method 
3 for calculating the slope. 

4. Parameter dependences 

For a fixed set of parameters which specify all the interactions, we study the depen- 
dence of 9r on the geometry, namely the aspect ratio (H/W) of the box and the linear 
size of the system. The parameters are kn = 10^, ks = 10^, 7„ = 5 x 10^, 7^ = 7„/100. For 
walls, kn is chosen to be 2 x 10^ to prevent particles to escape from the box. Since a sys- 
tem of particles with equal radii tends to form a hexagonal packing, we use particles with 
different sizes. The radii of the particles are drawn from a Gaussian distribution with the 
mean of 0.1 and width of 0.02, and the maximum (minimum) cut-off radius of a particle 
is 0.13 (0.07). In Fig. 2(a), we show the dependence of the angle Or on the height H, for 
values of /X = 0.2 and the width W = 2.0. Each angle is obtained by averaging over 20 
samples. The error plotted in Fig. 2(a) is the mean square sample-to-sample fluctuations. 
Here, we cannot see any systematic dependence on H. Also for other values of ji, we find 
that Or does not depends on the aspect ratio, as long as the ratio is sufficiently larger 
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than the slope of the pile generated. We then fix the aspect ratio to be 2, and study the 
dependence on the size of the system. The angle On for different values of W is shown in 
Fig. 2(b). All angles, as well as those presented in Fig. 2(b), are averaged over 20 samples, 
unless specified otherwise. For ji = 0.2, these angles decrease for small sizes, and seem to 
saturate starting around W = 3.5. For larger values of the angle saturates for larger 
values of W. For example, the angle continues to decrease until W — 4.0 for /j, = 0.3. On 
the other hand, for smaller fi, the angle saturates for smaller W. For jj, = 0.1, there is no 
obvious trend of the data, even up to the very small values ofW = 1.5. Since we want the 
angle obtained by this simulation not to suffer from a finite size effect, we will use in the 
foUowings the values W = A and H/W = 1 to calculate On. We will also study cases of fi 
not larger than 0.2. The simulation for larger values of /j, is limited due to the fact that 
one needs a larger aspect ratio and system sizes to be free of any finite size effects. 

Next, we study how On depends on the various interaction parameters in the system: 
7, k and fx. In a static configuration, the damping term is absent, so ideally 7 terms do not 
change Oji. However, since we prepared the sandpile by a dynamical method (by causing an 
avalanche), the angle may depend on 7. Also, kg is an the elastic constant of the artificial 
spring we introduced, so it should not make a difference in a static configuration, as long 
as we keep the value in a reasonable range. The quantity of particular interest is the 
friction coefficient since jj, determines if contacts between particles are stable ("stick") 
or unstable ("slip"). Since the stability of the whole structure (e.g. a pile) will be strongly 
influenced by that of individual contacts, we expect Or will be strongly dependent on fx. 
For example, 9r should be zero for = 0, if the individual grains in the pile are not 
moving. We first study the effect of kn and 7„ on Or. We will limit ourselves only to 
study the general trend such as the direction and the order of magnitude of the changes. 
We also fix // = 0.2. We measure the angle (inside parenthesis) for three different values 
of kn, 10^ (16.53 ± 0.73), 10^ (18.65 ± 0.70) and 10^ (17.99 ± 0.64). The difference in 
angle is very small, and there seems to be no systematic dependency. For three values of 
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7n = 50, 100, 500, the angles are 16.47 ± 0.75, 17.15 ± 0.62, 17.99 ± 0.64. The angle 
seems to decrease systematically, as 7„ is being decreased. However, the magnitude of the 
changes is still small (~ 5%). 

Now, we study the effect of /i. In Fig. 3, we show 9r obtained for several different 
values of In the range of /x we studied, there seems to be a linear relation between the 
angle and /x. This relation can be true for small values of /x, but it can not be true for the 
entire ranges of fi. The maximum 9ji is limited to 90 degrees, while the value of fi can 
be arbitrarily large. We will come back to discuss this relation later. Note that there are 
two friction coefficients in the system, one between particles and another between particles 
and walls. We will argue that, for a sufficiently large pile, the friction coefficient which 
determines du is that between the particles. Consider a sandpile on a table. The stress 
distribution near the top part of the pile would not be altered by the stress distribution 
at the bottom of the pile. Therefore, only the friction coefficient between particles can 
change Or in this region. On the other hand, the stress distribution near the bottom of 
the pile will greatly be influenced by the particle-wall friction coefficient. So, we expect 
the angle be different near the bottom of the pile, if the friction coefficients are different 
from each other. 

The piles discussed above are generated by causing avalanches. In experiments, the 
structure just after an avalanche (e.g. Fig. 1(c)) is not critical, but stable. In other words, 
one must apply an additional finite force to make the structure unstable. One way to apply 
the force is by tilting the box which contains the pile. The tilting angle 6t is defined as 
the rotation angle at which the pile becomes unstable. We want to emphasize that 9t is 
shown to be non-zero for real sandpile experiments. We measure the tilting angle for our 
model as follows. Starting from the pile like one in Fig. 1(c), we rotate clockwise the whole 
box with a constant rate of 10~^ degree/iteration. Then, we record the angle at which the 
pile starts to move, which is defined as the tilting angle 9t- The tilting angle for several 
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values of is shown in Fig. 4. Here, the width of the box W is 4, and the aspect ratio 
is 1. Indeed, one needs a finite tilting angle for the piles generated using our model, and 
it gives us confidence that the model studied here is reproducing the behavior of realistic 
static friction. The finite 9t implies that the pile is stable (not critical), therefore Oms 
should be larger than Or. 

5. Pile with constant flux 

In the previous section, we argued that On is smaller than 9ms based on the fact 
that the pile is stable even if it is tilted by a finite angle smaller than 6t- We now 
propose a method of obtaining the angle of marginal stability as well as the angle of 
repose. Consider an empty box without a right wall, and put one layer of particles at the 
bottom. We monitor the maximum velocity of the particles. If the maximum velocity is 
smaller than certain value Vcut, then we insert a new particle at the upper left corner of 
the pile. Once the particle is added, we then wait until the maximum velocity of particles 
is again smaller than Vcutj then add a new particle. This procedure is repeated long time 
for good statistics. In Fig. 5, we show the angle of the pile just before one inserts a new 
particle. Here, fi = 0.2, W = H = 4.0 and Vcut = 0.1. We also simulate the system with 
Vcut = 0.01, 0.001, and find no essential difference. The angle is zero at the beginning, and 
increasing until it seems to just fluctuate for iterations larger than 4 x 10^. The curve 
shown in the figure is quite noisy, which suggests that many configurations (or packings) 
are possible in the steady state. The maximum angle of the pile one can build up before 
avalanches is larger than the 9ft obtained before. This could be an additional evidence that 
our model reproduces the difference between 0ms cind 9ji. We can estimate the difference 
to be of the order of distance between two dotted lines in Fig. 5. Here, the dotted lines 
represents the mean square fluctuations of the angle. 



6. Theoretical approach 
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In the previous section, we measured various angles On, $t and $ms for the model 
sand. How can we understand these angles? In order to calculate these angles, one should 
know the stress field inside the pile. Liffman et af^ suggested an approximate way of 
calculating the field in a packing of equal sized spheres, which is illustrated in Fig. 6. In 
order to calculate the stress at a point O, one draws two lines of slope \/3 (line OA) and 
— \/3 (line OB) starting from the point O. Then, the length of the lines {I a and Ib) within 
the pile (and above the point) is approximately proportional to the force exerted by the 
pile. For a more detailed explanation as well as for the justification of this procedure, see 
Ref. 26. From this stress field, we calculate 9ms as follows. Consider a point at the bottom 
of the pile. The normal (to the bottom surface) force at that point is {Ia + is) sin(7r/3), 
while the tangential force is {Ia — Ib) cos(7r/3). We then apply the Coulomb criterion. If 
the ratio of the tangential to the normal force is larger than the contact is unstable. In 
this way, for given we obtain the range of angles at which the pile is stable. The largest 
angle at which the pile is stable is the angle of marginal stability, which is 

^MS = tan-i(3/x). (3) 

Since On + Ot is approximately the angle of marginal stability, we plot both measured 
Or + Ot and calculated 6ms by the above procedure in Fig. 7(a). The difference between 
the two angles is either due to the fact 9ms 9r + 9t or the error of the approximation. 
In fact, the approximation (and tt/S angle) is derived from an ordered packing of the 
particles with same radius. It is possible that the stress field in a disordered packing is 
very different from that of an ordered one. In that case, one needs a new approximation 
scheme to calculate the stress field. We also calculate 9t for given angle of repose 9ii and 
ji. It is given by 

_ i?cos(7r/3 - 9t) - cos(7r/3 + 9t) 
^ ~ i?sin(7r/3 - 9t) + sin(7r/3 + 9t) ' ^ 

where 

^ _ tan(7r/3 + 9t) + tan(gj^) cos(7r/3 + 9t) 
~ tan(7r/3 - 9t) - tan{9R) ' cos(7r/3 -9t)' 
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The 6t obtained by eq. (4) as well as the measured tilting angle are plotted in Fig. 7(b). 
One can also see that the difference between the two is small. Unlike the difference in 
Fig. 7(a), these two angles should coincide if the stress distribution is calculated correctly. 

The main point of presenting the theoretical approach is to show a "first" approximation for 
the problem, not to do quantitative comparison between the measured and calculated angle. 
In order to obtain more accurate numbers, one has to know a better way of calculating the 
stress field inside the pile. However, it is encouraging to see that even the values obtained 
by the first approximation are comparable to the measured ones. 
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Figure Captions 

Fig. 1: (a) Box filled with N = 1600 particles just before the right wall is removed. The 
thickness of the lines between centers are proportional to the normal force, (b) Average 
kinetic energy per particle (in erg) during the whole sequence of simulation. The 
energy initially increases as particles fall down, then decays with time. When the 
right wall is removed (iteration = 30000), it increases again, (c) Static pile obtained 
after the avalanche. 

Fig. 2: (a) The angle of repose 6ji vs the height H with the width W = 2 and /i = 0.2. The 
angle seems just fluctuate, and no systematic dependency is found, (b) The angle of 
repose 9^ vs W for several values of yu: = 0.1 (diamond), 0.2 (box) and 0.3 (circle). 
Here, the aspect ratio is kept to be 2. 

Fig. 3: The angle of repose Or vs /j, with W — H — A. 

Fig. 4: The tilting angle 9t vs /j, with W = H = 4. 

Fig. 5: The angle of pile measured when we add a new particle to the pile. The two dotted 
lines indicate the width of the fluctuation 

Fig. 6: The stress at a point inside the pile is approximately the vector sum of line A and B. 

Fig. 7: (a) The measured Or + 9t and the calculated 9ms is shown for several values of 
ji. There is a difference between the two. (b) The measured and calculated 9t for 
different values of The difference is smaller than that of (a). 
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